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Abstract: We construct a specific formalism for calculating the one-loop virtual cor- 
rections for standard model processes with an arbitrary number of external legs. The 
procedure explicitly separates the infrared and ultraviolet divergences analytically from 
the finite one-loop contributions, which can then be evaluated numerically using recursion 
relations. Using the formalism outlined in this paper, we are in position to construct the 
next-to-leading order corrections to a variety of multi-leg QCD processes such as multi-jet 
production and vector-boson(s) plus multi-jet production at hadron colliders. The final 
limiting factor on the number of particles will be the available computer power. 
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1. Introduction 



Multi-jet and vector-boson(s) plus multi-jet events are already common at the TEVATRON 
and will become even more frequent at the LHC. In addition to the information they 
contain about the strong dynamics, they also form important backgrounds to new physics. 
At present, the rates for such processes can be evaluated relatively easily at leading order 
(LO) in parton level perturbation theory [1, 2]. However, the predicted event rate depends 
strongly on the choice of the coupling constant (or equivalently the renormalisation scale) 
so that the calculated rate is only an "order of magnitude" estimate. In addition, there is 
a rather poor mismatch between the "single parton becomes a jet" approach used in LO 
perturbation theory and the complicated multi-hadron jet observed in experiment. While 
these problems cannot be entirely solved within perturbation theory, the situation can 
be ameliorated by calculating the strong next-to-leading order (NLO) corrections. This 
reduces the scale dependence and the additional parton radiated into the final state allows 
a better modelling of the inter- and intra-jet energy flow as well as identifying regions where 
large logarithms must be resummed. 

While many NLO calculations are available for 2^2 processes in hadron colliders 1 , 
relatively few estimates of 2 — > 3 processes based on five-point one-loop and six-point 
tree matrix elements exist. Currently, numerical programs capable of producing NLO 
predictions for fully differential distributions are available for 3-jet [4] and V+2-jet [5] final 
states, as well as the important pp — > tiH [6] and pp — ► ff+2-jet [7] signatures of the Higgs 
boson. The NLO predictions for important background processes such as pp — > yy+jet 
and pp — > QQ+jet are not available yet. 2 

The ingredients necessary for computing the NLO correction to a N particle process 
are well known. First one needs the tree-level contribution for the N + 1 particle process 
where an additional parton is radiated. Second, one needs the one-loop N particle matrix 
elements. Both terms are infrared (and usually ultraviolet) divergent and must be carefully 
combined to yield an infrared and ultraviolet finite NLO prediction. 

The radiative contribution is relatively well under control and can easily be auto- 
mated. Computer programs exist that can generate the matrix elements (and associated 
phase space) either Feynman diagram by Feynman diagram [2] or using recursion rela- 
tions [1] . The infrared singularities that occur when a parton is soft (or when two partons 
become collinear) can then be removed using well established (dimensional regularisation) 
techniques so that the "subtracted" matrix element is finite and can be evaluated in 4- 
dimensions [11]. 

In contrast, the one-loop contribution is usually obtained analytically on a case by case 
basis and extensive computer algebra is employed to extract the infrared and ultraviolet 
singularities in one-loop graphs that appear as poles in the dimensional regularisation 

1 Here, attempts are currently underway to evaluate the next-to-next-to-leading order (NNLO) correc- 
tions with the aim of reducing the theoretical uncertainty to the level of a few per cent - the anticipated 
experimental accuracy. See Ref. [3] for an overview. 

2 Some of the same five parton one- loop matrix elements [8] have been used to predict four jet event 
shapes in electron-positron annihilation [9] as well as 3 + 1-jet rate in electron- proton collisions [10]. 
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parameter e = (4 — D)/2 and then to evaluate the finite remainder in terms of logarithms 
and dilogarithms. At present, this is the bottleneck for producing NLO corrections for a 
variety of multiparticle processes. 

The NLO corrections to the 2^3 processes described above require one-loop cal- 
culations involving up to 5 external particles, i.e. the one-loop five-point function. The 
limiting factor is the number of external particles involved. While the one-loop A-point 
function can recursively be expressed in (A — l)-point functions all the way down to simple 
three-point functions [16], the algebraic complexity as N increases quickly becomes over- 
whelming. 3 This is mainly due to the rapidly increasing number of kinematic scales present 
in the problem. 

Adding one further particle to the final state corresponds to 2 — > 4 hadron collider 
processes. The NLO corrections to such processes requires the knowledge of one-loop six- 
point functions. To date, no standard model computations of one-loop six-point amplitudes 
exist. 4 

Rather than evaluating the one-loop matrix elements algebraically, an alternative 
would be a numerical approach [19, 20, 21]. There are two main problems. First, the 
ultraviolet and infrared singularities must be handled correctly so that they can be can- 
celled against those coming from the radiative contribution. Second, when anomalous 
thresholds occur in the loop integral, special attention must be given to the behaviour 
around (possibly) singular points in phase space [20]. 

Nevertheless, several avenues of attack are readily apparent. One approach that has 
already been successfully used for calculating five-point matrix elements would be to con- 
struct the loop amplitude by "sewing" tree amplitudes together [23]. Another way would 
be to combine the virtual and real contributions so that the divergences cancel directly in 
the integration over the loop momentum [24] . Alternatively, one can construct directly the 
counterterm diagram by diagram [25]. 

Here, we follow a different path. It is well known that (tensor) one-loop integrals can 
be written as combinations of finite four-point scalar integrals in D + 2 dimensions and 
infrared/ultraviolet divergent triangle graphs in D dimensions [12, 16, 22]. For example, in 
Ref. [26] the scalar six-point function is analytically expressed in terms of D dimensional 
triangles and D + 2 dimensional boxes using the recursive techniques of Ref. [16]. Using 
this result one can explicitly separate the divergent contributions, which are computed an- 
alytically in D dimensions using dimensional regularization, from the finite contributions 
which can be evaluated numerically in 4 dimensions. However, using the recursion rela- 
tions to calculate the divergent part of multi-leg loop integrals leads to an overwhelming 
algebraic complexity, particularly when A > 6 and when more kinematic scales are in- 
volved in the problem. While the infrared pole structure is always very simple, each level 

3 The fact that the one-loop TV-point function can recursively be expressed in (N — l)-point functions was 
well known for 4-dimensional integrals [12]. However, the extension to dimensionally regulated integrals 
was not obvious as the identities were based on the fact that space-time is four dimensional. The extension 
to higher dimensions was first formalized in Refs. [13, 14] and further developed by a variety of authors [15, 
16, 22]. 

4 Note that the six-point amplitudes have been evaluated in the Yukawa model [17, 18]. 
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of recursion introduces additional inverse powers of determinants of kinematic matrices 
which are difficult to cancel analytically. In this paper we advocate a different technique 
to directly extract the soft /collinear divergences from the loop graph and express them in 
the basis set of divergent D dimensional triangle graphs 5 . This approach to extract the 
divergent contribution allows us to minimize the algebraic manipulations needed to isolate 
the singularities. Once the coefficients of the divergent integrals are known, we can elimi- 
nate the divergent triangles from the basis set of integrals, leaving us only with a basis of 
(known) finite integrals. The kinematic coefficients multiplying the finite integrals can then 
be numerically evaluated using the recursion relations thereby maximizing the number of 
external legs that can be handled. 

Our paper is constructed as follows. First, in Sec. 2 we define the necessary formalism 
for relating tensor integrals in D ~ 4 dimensions to integrals in higher dimensions with 
additional powers of the propagators. We discuss the infrared and ultraviolet divergence 
conditions for these integrals and establish how to isolate the divergences using a set of 
augmented recursion relations based on the results of ref. [16]. The basis set that emerges 
is quite natural, divergent triangle functions and a collection of (finite) higher dimension 
box functions. We show how the recursion relations can be used to numerically evaluate 
the set of finite integrals. 

In Sec. 3, we show how to determine the singular contributions directly without use of 
the recursion relations. For simplicity we work in the limit where all of the internal particles 
are lightlike and construct a master formula for the divergent term using divergent triangles 
as the basis set of functions. The coefficients of each of the divergent three-point function 
is determined in Sec. 3.1.4. Our conclusions are summarized in Sec. 4. 

Several technical appendices have been added. Appendix A gives details about the 
derivation of the basic recursion relations together with explicit construction prescrip- 
tions for the numerical calculation of the kinematic coefficients appearing in the recursion 
relations. In Appendix B, we list analytic expressions for the basis sets of divergent in- 
tegrals that appear once the recursion relations have been applied are specified. Finally, 
Appendix C contains an explicit expression for the soft /collinear contributions from an 
on-shell massless A-point function. 

2. Outline of the Formalism 

In this section we give an explicit outline for calculating one-loop matrix elements with an 
arbitrary number of external legs. The strategy is to reduce the A-point one-loop integrals 
to a basis set of calculable master integrals. 

We start in Sec. 2.1 by relating the tensor integrals appearing in dimensionally reg- 
ulated Feynman diagram calculations (in D = 4 — 2e dimensions) to integrals in higher 
dimensions with additional powers of the propagators. This enlarged set of integrals nat- 
urally breaks down into three classes of scalar integrals: finite, ultraviolet (UV) divergent 
and infrared (IR) divergent integrals. We discuss the infrared and ultraviolet divergence 

5 Note that a similar approach has been employed by Dittmaier [27]. Although our method is very similar 
for the scalar integrals, it differs in the treatment of the tensor integrals. 
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conditions for one-loop integrals in arbitrary dimension and with arbitrary powers of the 
propagators in Sec. 2.2 and develop the necessary basis sets for the divergent integrals. 
In Sec. 2.3 we derive a complete system of augmented recursion relations needed for the 
formalism. Finally, in Sec. 2.4, we formulate the decomposition of the amplitude into a set 
of master integrals produced by the recursion relation. 

2.1 Decomposition of One-Loop Graphs to Master Integrals 

The one-loop M-particle amplitude Am with external momenta pi can be written as a sum 
over the individual one- loop Feynman graphs Mg, 

Am{Pi,P2,---,Pm) = ^2M g (pi,P2,...,Pm)- (2.1) 

G 

Each of the Feynman graphs is itself a sum over different rank-m iV-point functions, where 
m< N < M, 

N 

m g (pi,p2,...,pm) = EE c fiM 2 -^ 7 ^ 2 "' Mm ^;te},{i})- (2-2) 

G m=0 

The coefficient C^ lM2 ... Mm is typically composed of tree-level multi-particle currents that 
depend on the properties of the external particles (such as spin, polarization and ultimately 
on the momenta). In general this coefficient can depend on dimensional factors (i.e. it 
depends on the regulator e). Here, we are more concerned with the rank-m iV-point tensor 
integrals with unit propagators that are defined as 

T^Z^ST' (2 - 3) 

where the propagator terms are given by di = {I + qi) 2 + iO. 

As can be seen in Fig. 1, the momenta {%} that characterise the loop integral are 
composed from the external momenta. Depending on the topology of the diagram the qi 
are sums over different pi. The notation {1} indicates that all of the propagators are raised 
to unit powers. 

The generic tensor integral of Eq. (2.3) contains three different classes of divergences. 
The first class contains the so-called rescattering singularities. That is, two non-adjacent 
propagators are on-shell. These singularities are protected by the +i0 prescription and 
generate the complex part of the integral. In an analytic calculation these contributions 
are generated by analytic continuations of the transcendental functions into the physical 
region. The second class of singularities are generated when two adjacent propagators 
become singular. These are the genuine soft /collinear singularities. The final class of 
singularities occurs when the loop momentum becomes large - the ultraviolet singularities. 

To reduce the tensor integral to scalar integrals we use the method developed in 
Ref. [28] to extract the Lorentz structure and raise both the dimension and powers of 
the propagators within the remaining scalar loop integral such that 

I^-^(D;{ qi },{l}) = hx+E^i-m) xi\x 2 \---x N \ 

\,Xl,X 2 ,-.;X N 
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Figure 1: The generic TV-point loop graph. 



x /jvp + 2(m - A); {1 + a*}). 



(2.4) 



The structure {g^q^q^ 2 ' ' ' q^}^ 1 ^ 2 Mm means we assign the Lorentz indices in all distinct 
manners to a number of A metric tensors g a p, x\ momentum vectors q±, etc. For example, 



{qmr it12 

{gqiV 1 ^ 3 
{qlqlY 1 ^ 3 ^ 



(2.5) 



Finally, the generalised iV-point scalar integral with raised propagator powers is defined 
by 



I N (D- { Ql }, {v t }) = j 



d 



As we will see later, it is often convenient to classify the integral according to the value of 



a = v\ + v 2 H 1- v N . 



(2.7) 



Note that if one of the v \ = the corresponding propagator is removed yielding a (N — 1)- 
point function. Furthermore, for notational purposes we will often suppress the momentum 
arguments 6 



lN(D;{ qi },{^}) = I N (D;{^}). 



(2.8) 



6 We remind the reader that when a propagator is pinched out, and the topology of the loop integral is 
changed, the momenta {<?;} associated with that particular graph also changes. 
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Throughout the paper we will follow this notation which is similar to that developed in 
Ref. [22]. 

Generically, these scalar integrals in higher dimensions and with repeated propagators 
can be related using recursion relations, to three basis sets of integrals that reflect the 
singularity properties of the integrals: the set of finite integrals, J^ m , the set of infrared 
divergent integrals, X IR and the set of ultraviolet divergent integrals, T uv . The one-loop 
amplitude can therefore be written as, 

Am(j>i,P2, ■ ■ ■ ,Pm) = Affipi, . . . ,p M ) +A f M n (pi, . . . ,p M ) +Am / (pi, . . . ,p M ) 

= £ Kl R ll R {{ qj }) + ^ KM m ({ qj }) + £ Kf v lf v {{ qj )l 

it i 

(2.9) 

where the summations run over the integrals in each of the three basis sets. The kinematic 
factors K are functions of the kinematic scales in the process. To make Eq. (2.9) more 
concrete, we first classify each of the scalar integrals that can appear in Eq. (2.4) according 
to whether it is finite, IR divergent or UV divergent and then reduce it to the basis set 
using the appropriate recursion relations. 

2.2 Classification of Integrals 

The first important step is to identify the conditions under which a specific integral is finite, 
IR divergent or UV divergent. 7 Ultimately, these conditions will control which recursion 
relation to use in the reduction towards the master integrals. The integral In(D; {z^}) 
may be infrared divergent only when 

a-N + 2 = j, (2.10) 

depending on whether or not the external legs are on-shell and is ultraviolet divergent 
when 

D 

a = -. (2.11) 

a is defined in Eq. (2.7) and, since a > N, integrals may be finite even if the external legs 
are on-shell provided that D > 4 and the condition 

a-N + 2<^<a, (2.12) 

is satisfied. When an integral is IR and/or UV divergent, we wish to extract the divergences 
and relate the finite remainder to the set of finite integrals using recursion relations. 

First consider how the singularity conditions of Eqs. (2.10), (2.11) and (2.12) work 
when a and D vary for fixed N. This is illustrated in Fig. 2 for N = 2, 3, 4, 5. In each 
case, the IR divergent condition of Eq. (2.10) is shown as a red line, while Eq. (2.11) is 



7 In the following, we give a general discussion of the divergent conditions. A more detailed analysis of 
the singularity structure of integrals with iV < 6 including a thorough discussion of the role of internal 
masses and anomalous thresholds is given in Ref. [20]. 
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Figure 2: A plot in (D/2,a) for N = 2,3,4,5. The dots represent integrals that can occur in 
a Feynman diagram evaluation. The ultraviolet divergent integrals lie along the lower solid line, 
while the infrared divergences lie along the upper solid line. For N = 2, the IR and UV divergences 
coincide and are shown as a magenta line. Finite integrals lie on the dashed lines. The ultraviolet 
divergent six-dimensional bubble graph (produced by the second rank tensor bubble) is denoted by 
a blue circle-cross. 



shown as a blue line. For bubble graphs (Fig. 2(a)), the IR and UV lines coincide and are 
drawn in magenta. The dots represent the integrals that appear when Eq. (2.4) is applied 
to tensor integrals up to rank N. In all cases, the IR line passes through origin, a = N 
and D = 4, indicating that the D = 4 scalar integral may be IR divergent. We also see 
that the allowed triangle integrals are always either IR divergent (unless the legs are all 
off-shell) or UV divergent. On the other hand, for N > 4, families of finite integrals appear 
starting with the scalar integral in D = 6. 

It is also instructive to consider the finiteness of integrals of varying N with fixed 
dimension. This is shown in Fig. 3 where the IR and UV conditions are again shown red 
and blue while the dots correspond to integrals produced by Eq. (2.4). Fig. 3(a) shows that 
all four-dimensional integrals can be divergent while the first integrals that are guaranteed 
to be finite independent of whether the legs are on-shell or not are higher dimension box 
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D=4 



D=6 




Figure 3: A plot in (N,cr) for D = 4,6,8, 10. The dots represent integrals that can occur in a 
Feynman diagram evaluation. The ultraviolet divergent integrals lie along the blue horizontal line, 
while the infrared divergences lie along the red line. Finite integrals lie on the dashed lines. The 
ultraviolet divergent six-dimensional bubble graph (produced by the second rank tensor bubble) is 
denoted by a blue circle-cross. 



graphs. 

Reading off the allowed UV divergent integrals from Fig. 3, we see that a suitable (but 
overcomplete) basis set is 

l uv = {I 4 (D = 8 - 2e; 1, 1, 1, 1), h(D = 8 - 2e; 2, 1, 1), I 2 (D = 8 - 2e; 3, 1), 
I 2 (D = 8 - 2e;2, 2),I 3 (D = 6 - 2e; 1, 1, 1),I 2 (D = 6 - 2e; 2, 1), 
h(D = 6 - 2e; 1, 1), I 2 (D = 4 - 2e; 1, 1)} (2.13) 

such that 

A M (PUP2, ■ ■ ■ ,PM) = £ Kl R ll R {{q 3 }) + £ Ki^dgj}) + E ^r*F y ({<?,-})■ 

(2.14) 

While the UV divergent functions are all known analytically [14, 29, 30], the basis sets of 
finite and IR divergent integrals have an infinite number of members. To derive finite basis 
sets of known integrals for the latter two classes we employ recursion relations. 
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2.3 Recursion Relations 

In this subsection we show how to map the higher dimensional finite and IR divergent 
integrals onto a basis set of known finite and divergent integrals. This will fully specify 
Eq. (2.9) and produce the explicit set of master integrals {2 UV ,l IR ,lf m }. 

Because the external momenta are 4-dimensional, the recursion relations come in two 
distinct classes that depend on the number of external legs. When the number of legs is 
five or less, the external momenta form a vector basis for the resulting tensor structure. 
However, for more than five external lines, the set of external vectors is over-complete and 
the net result is that the recursion relations change for N > 6. 

First, we introduce some notation and define the symmetric kinematic matrix 

Sij = (ft ~ Qj) 2 - (2.15) 

By adding a row (or column) over the inverse of the kinematic matrix we make the further 
definitions (provided the inverse of the kinematic matrix exists), 

3 i ij 

The Gram matrix 

Gij = 2 qi • Qj , (2.17) 

is closely related to the kinematic matrix. 

The recursion relations needed to relate all possible integrals to the basis sets are for- 
mulated as follows: 

I. N < 5: det(S) + 0, det(G) ^ 

We will use four recursion relations in this case, one basic recursion relation and three 
derived, composite recursion relations. 

The basic recursion relation is obtained by integration by parts as explained in ap- 
pendix A. It is given by, 

[y k - l)I N (D;{ui}) 

N 

= ~ E S k^N(D - 2; fa - S H - S lk }) -b k (D- a) I N (D; fa - 5 lk }), 

i=l 

(2.18) 

where a is defined in Eq. (2.7). This relation either reduces both the dimension D and the 
accompanying value of a by two, or it keeps D fixed and reduces a by one unit. Its action 
is indicated by the red lines in Fig. 4. This identity is vital for extracting the infrared 
singularities from the integrals produced by the Davydychev decomposition of Eq. (2.4). 

Note that for the case iV = 3 in Fig. 4(a) we assume that det(S) ^ which is only 
the case if all three external momenta are off-shell. If any of the external momenta is an 
on-shell lightcone vector, the inverse of the kinematic matrix does not exist and Eq. (2.18) 
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is invalid. For practical purposes, this is not a problem because analytic expressions for 
arbitrary dimension and powers of propagators exist for triangle integrals with at least one 
on-shell leg and massless internal propagators (see Appendix B). In four dimensions, the 
triangle integral with three off-shell legs and unit propagators is finite and this serves as a 
reminder that integrals on the "IR" lines in Fig. 4 are not necessarily IR singular, rather 
they can be IR singular if some of the external momenta are lightlike. 

Eq. (2.18) is a more detailed equation than usually quoted in the literature and it can 
only be applied to scalar integrals with at least one index v k > 1. However, we can use it 
to derive the standard recursion relation [16], 

N 

{D — I — a) B I N (D; fa}) = I N (D - 2; fa}) - £ bJ N (D - 2; fa - 5 H }). (2.19) 

i=i 

Details of the derivation can be found in Appendix A. Eq. (2.19) reduces D by 2 and a by 
1. Its action is illustrated by the blue lines shown in Fig. 4. For the case that all vi = 1 
and D = 4 — 2e, we will use this recursion relation to extract the infrared singularities from 
graphs with at least one on-shell leg as pinched integrals together with scalar integrals in 
D = 6 — 2e dimensions, 

N 

In(D; fa}) = (D + l-a) B I N (D + 2; fa}) + £ hI N (D; fa - 5 U }). (2.20) 

i=l 

Its action is illustrated by the magenta line in Fig. 4. 

By combining Eqs. (2.18) and (2.19) we can eliminate the dimensional prefactor to 
yield the fourth recursion relation, 

fafc-lJJjvGDjfa}) 

= - 2; fa - S lk }) + jz( bj §- S m) MD - 2; fa - S H - %}). 

i=i ^ ' 

(2.21) 

This identity reduces D by 2 and a by either 1 or 2 units and its action is indicated in 
Fig. 4 by the green lines. 

When N = 5 it is easy to avoid the UV line because of the presence of two types of 
finite integrals (D/2 = a — 1 and D/2 = a — 2). However, this is not the case for N = 3 
and N = 4 and we are forced to use Eq. (2.21) which appears to produce UV divergent 
integrals from finite integrals. This implies that in the sum over the produced UV divergent 
integrals the divergences must cancel. Indeed, because 

i=i v 7 

these artificial UV divergences in the individual integrals cancel in the sum as expected. 
This means we replace each UV divergent integral in Eq. (2.21) by a UV finite regulated 
integral, 

F N V {D = 2a; fa}) = I U N V (D = 2a; fa}) - R(D, a), (2.22) 
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Figure 4: A plot in (D/2,a) for (a) the finite triangle graph, N = 3, (b) the box graph, N = 4 
and (c) the pentagon graph, N = 5. The ultraviolet divergent integrals lie along the lower solid 
line, while for N > 3 any infrared divergences lie along the upper solid line. Finite integrals lie on 
the dashed lines. The action of Eqs. (2.18), (2.19) (2.20) and (2.21) are shown red, blue, magenta 
and green respectively. 



where the scaleless regulator term depends on the dimension and a with the constraint 
D = 2a, 

After the subtraction the integral is rendered finite. Of course, for this replacement to be 
useful we must require that Eq. (2.19) is still valid for the regulated integrals In. Therefore, 
the regulator function must obey 

JV 



(D - 1 - a) B R(D, a ) = ~Y^ hR(D - 2, a), (2.24) 
or equivalently 

(D-cr- l)R(D,a) = -R(D-2,a), (2.25) 



which is trivially verified. 

The application of the various recursion relations is clear from the figures: 

• If N ^ 3 and D = 4 apply Eq. (2.20) (magenta arrow) 
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• If N ^ 3 and integral is on the IR line (D/2 = a - N + 2) apply Eq. (2.18) (red 
arrows) 

• In the other cases one can apply Eqs. (2.19) (blue arrows) and (2.21) (green arrows) 
with the restriction a — N + 2 < D/2 < a. Any "fake" UV divergent integrals 
generated can be regulated as explained in the previous paragraph. 

These simple rules are sufficient to reduce all integrals with N < 5 and non-exceptional 
kinematics to the basis set of integrals. 

II. N = 6: det(S) + 0, det(G) = 

For N = 6 the kinematic matrix Sij is still invertible. This means that the same 
recursion relations hold as in the A < 6 case. However, because det(G) = 0, B = J2i hi = 
[16]. This does not affect the basic recursion relation given in Eq. (2.18), but does alter 
the composite recursion relations and Eqs. (2.19)-(2.21) all collapse to yield 

6 

h(D- {ui}) = biWD; {n - 5 H }) . (2.26) 
i=i 

This equation preserves the value of D, but reduces a by unity. Its action is indicated by 
the green arrow in Fig. 5. An additional equation that reduces D is also required and can 
be obtained by combining Eqs. (2.18) and (2.26), 

b k (D-a + v k -l)I 6 {D;{vi}) = 

6 6 

-u k hh(D; {vi + S lk - So}) - S klh{D - 2; {v x - 5 H }). (2.27) 

i^k i=l 

This recursion relation works in 2 different manners. The second term reduces D by 2 and 
a by unity and acts in a similar manner to the identities for N < 5 case and is illustrated 
by the magenta arrow in Fig. 5. On the other hand, the first term preserves both D and 
a while reducing N. 

The application of the recursion relations for this case are straightforward. 

• If TV = a apply Eq. (2.26). 

• If D/2 = a - N + 2 (i.e. on the IR line) apply Eq. (2.18) (red arrow). 

• If D/2 < a — 1 apply Eq. (2.26) or (2.18) (green or red arrow). 

• If D/2 = a — 1 apply Eq. (2.27) (magenta line). The index k is chosen such that 
v k = max;(z^) (this is not necessarily a unique choice). 

These simple rules are sufficient to reduce all integrals with = 6 and non-exceptional 
kinematics to N = 5 integrals. 
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Figure 5: A plot in (D/2,a) for hexagon graphs, N — 6. The ultraviolet divergent integrals 
lie along the lower solid line, while the infrared divergences lie along the upper solid line. Finite 
integrals lie on the dashed lines. The action of Eqs. (2.18), (2.26) and (2.27) are shown red, green 
and magenta respectively. Note that integrals produced by pinching so that N < 6 are not shown. 

III. N > 7: det(S) = 0, det(G) = 

In this case, the inverse of the kinematic matrix does not exist. This invalidates the 
use of Eqs. (2.18)-(2.21) for N >7. 

However, the fact that the kinematic matrix is made up of 4-dimensional vectors which 
obey momentum conservation leads to additional relations. To derive these new types of 
recursion relations valid for these high values of N, we use a generalized inverse based on 
the singular value decomposition of the Gram matrix. 

The first basic recursion relation is based on the kernel of S. By solving (S.z)i = for 
z we obtain [22] 

z k I N (D; fa}) = -J2 ZiIff(D; fa + 5 lk - 5 U }), (2.28) 

or alternatively, 

zJ N (D; fa + % - S u }) = 0, (2.29) 

i 

where the explicit Zi are given in Appendix A. An important property is that Y2i z i = 0- 
Eq. (2.28) keeps D and a constant, but can pinch out particular propagators to reduce N. 
Its effect is illustrated by the magenta arrows in Fig. 6. 

The existence of the second recursion relation is highly non-trivial. It is based on the 
fact that the unit vector Jj = (1,1,..., 1) lies in the range of the singular matrix S, i.e. 
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D=16 




UV 



N 



Figure 6: A plot in (N, a) for D = 16 illustrating the procedure for evaluating integrals with 
N > 6. The dots represent the integrals produced in the Davydychev decomposition of rank 6, 7, 8 
and 9 tensor integrals. The ultraviolet divergent integrals lie along the lower solid blue line, while 
the infrared divergences lie along the upper solid red line. The action of Eq. (2.28) and Eq. (2.30) 
are shown as magenta and cyan respectively. 



the equation (S.r)i = U has a solution for rj. The existence of the solution and the explicit 
construction of the vector n are detailed in Appendix A and leads to [22] 

I N (D; fa}) = nI N {D; fa - 5 H }) (2.30) 

i 

where Y^i r i = 0- While in functional form this recursion relation is equivalent to Eq. (2.26), 
its origin is rather different. This is reflected in the fact that the unique coefficients bi = 
Y^j S~j l are unrelated to the many possible and equivalent sets of r, coefficients. For 
example, replacing — > + az; L leaves Eq. (2.30) unaltered. Eq. (2.30) also preserves D, 
but reduces a (and possibly N) by unity as illustrated by the cyan lines in Fig. 6. 

The application of the recursion relations for this case are similar to the N = 6 case: 

• If N = a use Eq. (2.30) 

• If D/2 < a — 1 apply Eq. (2.30) (cyan arrows). 

• If D/2 = a — 1 apply Eq. (2.28) (magenta line). The index k is chosen such that 

= m&xifa) (this is not necessarily a unique choice). 

These simple rules are sufficient to reduce all integrals with N > 7 and non-exceptional 
kinematics to = 6 integrals. 
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2.4 The basis set of integrals 



The net result of applying the recursion relations is that any amplitude can be written as, 

A M (PUP2, ...,PM)= Yl K v*2vA R ( D = 2 (° ~ I)! ^ ^2, V3) 

+ J2 Kf { M i N V (D = ^{-e}) 
+ Yl K[ m l( m (D = 4; 1,1,1) 

triangles 

+ ^K( m l( m (D = 6; 1,1, 1,1) 

boxes 

+ K[ m l[ m (D = 6; 1,1, 1,1,1) 

pentagons 

+ X>r^r- (2.3i) 

i=i 

Note that the absence of new basis integrals for N > 6 is intimately related to the dimen- 
sionality of space-time (i.e. det(G) = for N > 6). All of the integrals appearing in this 
basis set are analytically known. The IR and UV divergent integrals with massless internal 
propagators are known analytically for arbitrary powers of the propagators and arbitrary 
dimensions [29, 30]. For the sake of completeness, we list them in the Appendix B. Note 
that the regulated UV divergent triangles are produced by the action of Eqs. (2.21) 
and (2.26) and the sum is therefore finite. 

As can be seen in Fig. 2, UV divergent integrals occur for 2-point, 3-point and 4-point 
functions. By applying Eq. (2.19) to UV divergent 4-point functions, the UV divergences 
are moved to 3-point functions. This means the UV basis set of Eq. (2.13) is reduced to 
only 2-point and 3-point functions. Analytic formulae for these UV divergent integrals are 
listed in Appendix B. 

The finite 4-dimensional triangle integral with off-shell legs, l[ in , and the 6-dimensional 
box graphs are also known in terms of logarithms and dilogarithms [14, 29]. Furthermore, 
it is an empirical fact that in final expressions for physical quantities, the 6-dimensional 
pentagons do not appear [14, 16, 22]. This means K[ m = 0(e) and we do not need to 
know h(D = 6; 1, 1, 1, 1, 1) for NLO calculations. 

3. Determining the Infrared Divergent Contributions 

In this section we derive the techniques for determining the coefficients K IR that multiply 
IR divergent triangles l£ R in Eq. (2.31). All of the other coefficients that multiply the finite 
integrals can be determined numerically using the recursion relations. However, because the 
coefficient K IR multiplies a divergent integral, we need to know it in D dimensions. This 
can only be achieved using analytic methods. In Sec. 3.1, we decompose the divergent part 
of a rank-m iV-point function into a sum of tensor structures multiplied by the divergent 
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triangles and a kinematic factor, C N ' J . This kinematic factor is the building block for the 
coefficient K IR . 

First we derive an analytic expression for C 1 ^^ in Sec. 3.1.4. As it turns out this 
coefficient does not depend on the rank of its parent integral and all the tensor information 
is carried by a kinematic Lorentz structures. 

While in this section we consider only massless internal lines, the generalization to 
include masses is straightforward and has been outlined in Ref. [27]. By combining results 
from [27] with the techniques to deal with tensor integrals explained in this section, one 
readily obtains the extension to the massive case. 

3.1 Derivation of the Divergent Coefficients 

As we saw in Sec. 2, the IR divergences of any (tensor) loop integral can be isolated as a 
combination of triangle integrals with one and two light-like legs [16, 31]. In the standard 
approach using the recursion relations, the IR divergent graphs are produced from many 
different sources and ultimately drop out at the end of a lengthy algebraic calculation, 
often after large cancellations between terms. Here it is our aim to calculate the kinematic 
coefficients multiplying these infrared divergent integrals ab initio without recourse to the 
recursion relations. In other words, starting with some given (tensor) loop integral, we 
would like to take the IR limit directly and isolate the singularities, 

I^-^iD; vi,... v N ) IR -^? it Sft'"*", (3.1) 

where Sn contains all the infrared singularities. Ultimately, this will produce exactly the 
correct coefficients Klft 2 „ 3 that would be produced by application of the recursion relations 
so that 

N 

G m=0 

= E K™ 2U3 4 R (D = 2(a - 1); v u U2, u 3 ). (3.2) 

As we will show, identifying S^ '^" 1 is straightforward for scalar integrals but is slightly 
more involved for tensor integrals. 

3.1.1 Scalar Integrals 

The infrared singularity corresponds to the limit where the incoming loop momentum at 
the vertex where particle i joins the graph becomes collinear with the incoming momentum 
Pi. If we define t- L = t + p\ + ■ ■ ■ + Pi-i (&\ = I) to be the momentum entering the vertex 
and ii + pi to be the momentum exiting the vertex, then the collinear limit occurs when 
ti — > xpi. Both propagators adjacent to the incoming momentum pi are now on-shell, 
£ 2 = x 2 p 2 = and (£i + Pi) 2 — > (1 + x) 2 pf = 0, and produce infrared singularities. We 
therefore can construct a term that matches all of the infrared singularities when the loop 
momentum becomes collinear with one of the external momenta. With the identification 
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dN = do the counterterm for a scalar TV-point function in the limit £2 —* %P2 is given by 



S N ,2 = A( 



m D l 2 d\d2 \_dz{x)d/t(x) ■ ■ ■ dx{x) J ^ 



(3.3) 



*xp 2 



where A(p^) = 1 when p\ = and zero otherwise (i.e. when P2 is off-shell). The square 
bracket is to be evaluated in the limit that I2 — > xp2, or equivalently £ — > xtft — (1 + 
So far, the value of x has not been determined. However, there are still singularities when 
x takes on the specific value X — Xj such that dj(xj) = 0. These singularities can be 
extracted one at a time so that, 



N 



S N , 2 = A(pl)£l> 12J (1,1,1) 4 2j , 
i=3 



where, 



P 12J (1,1,1) =I 3 (£ = 4-2e; {q u q 2 , qj}, 1, 1, 1) = ^ 



m D l 2 did 2 dj 

is independent of the number of legs in the graph. The kinematic factor is given by, 

1 



N 



n^3 di(xj) 



(3.4) 



(3.5) 



(3.6) 



We now have accomplished our aim for the scalar integral. The IR singularities associated 
with incoming light-like momentum P2 is a sum of kinematic factors multiplying scalar 
three-point master integrals. 

3.1.2 Tensor Integrals 

In the case of tensor integrals the derivation is a bit more complicated. In the collinear 
limit £2 — > xp2 we find 



S%"-»» = A(p 2 2 ) 



dPi £^ m 3 



ITT 



D/2 



did 2 



£fJ-m 3 + l . . . g^k 

d^{x)d^{x) ■ ■ ■ djy(x) 



(3.7) 



IR limit 



where of the m loop momenta are integrated over and the remaining (m — 777,3) re- 
evaluated in the IR limit £2 — > xp 2 - A priori, any value of 7773 between and m is possible 
and leads to the correct IR behaviour. In Ref. [27] it is proposed that 1713 is chosen to 
be equal to m. Using the recursion relations of Sec. 2 this particular choice would give a 
rather complicated structure of IR triangles. However, as explained above more generic 
choices for the tensor structure or the IR terms are possible. 

To decide on the appropriate partition of loop momenta, it is important to consider 
the structure of the IR triangles produced by the recursion relations of Sec. 2.3. Once we 
can get the correct correspondence between the two different methods for evaluating the IR 
contributions we have by principle of consistency found the right procedure for partitioning 
the loop momenta. 
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13 1 




D/2 

Figure 7: A plot in {D/2, a) showing how the IR divergences from high-rank iV-point integrals are 
mapped onto IR divergent triangle integrals using Eqs. (2.18) and (2.30). The solid lines show the 
IR trajectories for TV = 3, . . . , 7. The points labelled (TV, m) correspond to integrals produced by 
rank-m TV-point tensor integrals. The red arrows show how Eqs. (2.18) and (2.30) systematically 
reduce both TV and a until the integral ends on the TV — 3 IR line. 

Fig. 7 shows how the IR divergences from rank-m TV-point integrals are mapped onto 
IR divergent triangle integrals using Eqs. (2.18) and (2.30). Because the recursion relations 
change in a fundamental matter once TV > 7 we have to distinguish these two cases. 

For TV < 6, systematic application of identity (2.18) reduce both a and D in such a way 
that the integral jumps from the TV-point IR trajectory to the (TV — l)-point IR trajectory. 
Each jump reduces a by 2 and D/2 by unity. Repeated application of the identity will 
always end on the N = 3 IR line. Provided that the rank-m satisfies m < V — 3, then we 
obtain triangle integrals with a = 3. However, if m > V — 3, then the recursion relations 
hit the IR triangle line at a > 3 and D/2 > 2. In Fig. 7 this is illustrated for rank-5 
six-point integrals (6,5). After reduction, we are left with triangle integrals in D/2 = 4 
and with a = 5, precisely the integrals that are produced by rank-2 triangle graphs in 
D/2 = 2. 

For N > 7 the recursion relation (2.30) works rather differently. It keeps D/2 fixed, 
but reduces a and V. Repeated application will eventually move the integral onto the 
V = 6 IR line and the previous case is applied. 

This leads us to the surprisingly simple conclusion that the correct partitioning of loop 
momenta in Eq. (3.7). For TV < 6 we find 7713 = max(m + 3 — A^, 0), while for TV > 7 we 
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have 777,3 = max(777 — 3, 0). Or summarized in one formula 



7773 = max(777 + 3 — min(iV, 6), 0). (3.8) 

When N < 6 the maximum value of 7773 is only three (reached for a maximum rank-777 = N 
iV-point function). However, if N > 6, the maximum value of 7773 is N — 3. Note that only 
the IR divergent part of the tensor triangle integral is required. 

Inserting this value of 7773 in Eq. (3.7) and partitioning the loop momenta in a sym- 
metric way the counterterm can be written as 

N 

S%y™ = A(pi) C^T(rr" - , (3.9) 

j=3 



where 



The notation 



rr-fll—flm T)^ 1 " '^ m 3 -pA t m 3 +1 ' • -fJ-m] /r> -. p.\ 

1 12,j - U 12,j ' 12j • ^■ LU) 



^ k\(m- k)\ 

shows how the Lorentz indices are distributed between the terms and includes a sum over 
all permutations of the indices. The kinematic factor is defined in Eq. (3.6) while the 
tensor structure remaining in the IR limit is given by, 

r£p = [ti...t*}^ XjP2 _ pi . (3.i2) 

Finally the IR part of the rank-/c triangle function is given by, 

_ ( n Xl„X 2 k-X 1 -X 2 \ Vl "' Vk 

X1X2 

x4 R (D = 4 + k - 2e; q u q 2 , qj, m + 1, x 2 + 1, k + 1 - xi - x 2 ). (3.13) 

Eq. (3.9) produces the IR singularities for the rank-m iV-point integral in the l 2 — > xp 2 
limit. Note that only those 3-point functions which are divergent (i.e. have at least one 
massless on-shell leg) are kept. 

We see that Eq. (3.9) contains exactly the same IR triangles as would be obtained by 
using the formalism of Sec. 2 to a rank-rn iV-point function followed by application of the 
recursion relations. Therefore, by inspection, the divergences of the Eq. (3.9) are identical 
to the original tensor integral. Because the kinematic factor C^' J does not contain any 
dimensional term (i.e. does not depend on the regulator e), the finite parts will also match 
exactly and, 

" = 0. (3.14) 



lim 

h— >XP2 



TfJ.l---fJ.rn _ qfJl---fJn 
A N °N,2 
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3.1.3 Sum over all IR limits 



We have constructed the IR divergent part of the loop integral associated with the partic- 
ular limit £2 — ► xp2- Similar IR singular terms for the other cases (i.e. ii — ► xpi) can be 
obtained by cyclic permutations of Eq. (3.4). That is (£i,pi,di) — > (^2,^2,^2) —>•■■—> 
{£n-iPn, cZjv) — > denoted by C(12---n). However, care has to be taken when 

adding up all contributions to construct the overall IR singular term. As is well known 
for real emission soft/collinear divergences, adding up all collinear limits will double count 
the soft double singularities. The same occurs for the virtual soft poles when two adjacent 
external momenta are on-shell. For example, when both P2 and p 3 are on-shell we get an 



identical soft/double pole contribution from Sn^ and Sn'- These terms are respectively 

12 3 23 1 

C n ' ld\did^ and C n ' /(fad^di. To correct for the double counting we have to subtract this 
singularity. With this in mind we can now construct the overall IR singular term for the 
loop integral, 



:(3) 



- E [• 

C(12-n) 

= E 

C(12-n) 

= E 

C(12-JV) 



°JV,2 



A(pi)A(pi)T^ C 



,12,3 

N 



3=3 



2\/T-All-"Mr 



A(pi)A(pi)Tf 2 J^ C 



N 



7l2,j 
"AT 



2 a(pI)^;-^ c 
+A(pl)A{pl)rr 2 K^e 



12,3 

N 



+A(p|) (l-A(p§))7S^C 



12,3 

N 



+A(pl)(i-A(pl))T£; 1 -^c 



23,1 
TV 



(3.15) 



Eq. (3.15) clearly shows the different contributions. Just as for real emission soft/collinear 
behavior, there is a double pole soft contribution when two neighboring incoming momenta 
are on-shell, A(p\) = A(pf) = 1. If one of the two neighboring incoming lines is off-shell 
there is a single pole collinear contribution. However, unlike the real emission case, there 
is an additional single pole contribution from non-adjacent propagators. 



3.1.4 Evaluation of C 



12j 
N 



If the momentum P2 is light-like, then the condition that £2 — > xp2 means that all of the 
propagator factors become first order polynomials in x. In this case, we can use, 



?2+P2 + --- +PjY 



(1 + x)s 2 -j - xs 3 ...j,= Aj(x - Xj 



(3.16) 



with 



A, 



82-j ~ S 3 ...j, 



(3.17) 
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x = s Jzi (3.18) 

S 3 ...j - 8 2 ...j 

Sk-j = (Pk H ^Pj) 2 5 ( 3 - 19 ) 

so that 

JV 



j=3 



The coefficient is merely the residue evaluated at the singular point x = Xj , 

Cn J = — 1 • ( 3 - 21 ) 

Wm=3,m^j A m (Xj — X m ) 

Inserting the definitions of Ai and Xi given in Eqs. (3.17) and (3.18) we find 

I \N-3 

r u,j y- S2 -j z mtil ci o<?\ 

llk=3,k^j ( s 2-jS3-k ~ S 3 :.jS 2 ...k) 

where % = YTj=iPi such that the formula respects cyclic permutations of the indices. 
Eq. (3.22) holds whether or not p 2 = 0. 

The special cases associated with neighboring incoming momenta are obtained by 
setting j = 3 in Eq. (3.22) and noting that S3. ..3 = p\ so that, 

M2,3 _ ( s 23 ~ P3) ,„ 9 ox 

llfc=4 («23S3-fc -P3 S 2-k) 

corresponding to evaluating t — > X3P2 — qi with X3 = S23/(-S23 — P3) and q\ = p\. The 
companion coefficient is given by making the cyclic momentum permutation, p\ —> p 2 ,P2 — * 
P3, Qi ->■ Q2 on Eq. (3.22) so that s 2 ...j -» S3...1 = pi, s 3 ...j -> S4...1 = S23, 



,23,1 _ ( g 23 ~ Pi) 

11^=4 {s23S3-k ~ plsA-k) 



Cf l = ^ V ^ FV , „ (3-24) 



corresponding to £ — > X1P3 — pi — P2 with xi = — p 2 / (p 2 — S23) 
When p 2 = P3 = the two coefficients become identical 



,,12,3 
AT 



^23,1 



(3.25) 



As an application of Eq. (3.15), in Appendix C we write down an expression for the 
IR singularities for a scalar massless A-point function with all external momenta on-shell. 



4. Outlook and Conclusions 

The formalism outlined in this paper offers a realistic prospect of constructing NLO Monte 
Carlo programs for processes involving high multiplicity final state particles. The key result 
is a strategy for combining the analytic evaluation of the IR and UV singular contributions 
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to loop integrals by using a suitable basis of divergent integrals with a system of recur- 
sion relations which allows the numerical determination of the finite integrals. This mixed 
approach avoids the usual algebraic log-jam in loop calculations that occurs when many 
particles and kinematic scales are present. In this paper, we have focussed on massless 
particles circulating in the loop. However, the extension to massive particles is straight- 
forward and opens up the possibility of computing the NLO corrections to any Standard 
Model process. 

The next phase of this project is the actual numerical implementation of the algorithms 
such that a rank-m iV-point function can be calculated. There are two avenues of approach. 
The first approach is to use the recursive algorithms numerically. The advantages are 
simplicity and all kinematic coefficients are calculated numerically. This allows us to check 
that e.g K 5 = and all K IR can be compared to the analytic calculation of these divergent 
coefficients. The recursive algorithm is a good diagnostic tool. 

However, if the aim is to calculate a full amplitude made up from many graphs con- 
taining different integrals of varying rank, the recursive approach would not be optimal. 
This is because the calculation of a rank-m iV-point function uses some of the calculations 
needed for evaluating a rank-/c M-point function (k < m, M < N). To avoid duplicate 
calculations and irrelevant paths, one could use a constructive method for evaluating the 
integrals. The algorithm starts by calculating all the relevant finite basis integrals. Then, 
step by step, we calculate the integrals with higher D, N and a using the recursion rela- 
tions in the reverse sense. For example, assume all Im{D;{vi}) are known with a < sq 
(i.e. each integral is a number in a look-up table). Using the recursion relations of Sec. 2, 
we can then calculate all the In(D; {yi}) with a = sq + 1. Now we know the look-up table 
for all integrals with a < sq + 1. We continue with this procedure until all integrals with 
M < N are known. After that, building the amplitudes is easy. All the tensor structures 
from the Davydychev decomposition are contracted in with the tree level currents. We can 
simply look up the numerical value for the appropriate integral In{D; {v{\) in the table. 
This method is very efficient. However, some of the diagnostic tools have been lost. (E.g. 
we assume K§ = and all divergent coefficients are correct). 

Taken together, the combination of recursive numerical and analytic algorithms should 
provide a compact and efficient way of evaluating multiparticle one-loop amplitudes, thereby 
opening up the possibility of estimating the NLO QCD corrections for processes such 
as pp — > N jets in association with vector-boson(s) or Higgs particles. Or for example, 
pp — > QQ plus N jets at NLO, pp — > QiQi + Q2Q2 plus N jets at NLO, etc. The only 
limiting factor on the particles is that the spin/helicity is less than or equal to one and 
(numerical) computing power. We expect that this latter problem will be solved with the 
advent of the LHC Computing Grid. 
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A. Derivation of Recursion Relations 



The one-loop recursion relations are based on the integration by parts techniques developed 
in Ref. [33] for calculating higher order beta functions. The application to one- loop multi- 
leg recursion relations was developed in Ref. [12, 22]. Here we will initially follow Ref. [22] 
by applying the integration by parts identity to 



(A.l) 



N 



After some trivial algebra and applying the dimensional shift identity 

N 

I N (D - 2; {v k }) = -J2»i MD; W + S ik }) , (A.2) 
i=l 

one finds the base equation 

N / N \ 

Y Y S i iVi u J i n(D; {vi + <%}) 
j=l \i=l ) 

N / N \ / N \ 

= - VilN(D - 2; {v x - S H }) - ID - 1 - "j Vi) W ' ^» ' ( A - 3 ) 

i=l \ j=l J \i=l J 

which is valid for any choice of the N parameters {yi}. 

To reduce tensor integrals in a systematic way, more general recursion relations are 
needed than those previously given in the literature [12, 22]. For any useful identity, the 
parameters {yi} must be chosen such that the matrix expression Yli=i SjiUi is either a delta 
function, zero or a unit vector. In all these cases, the left hand side of Eq. (A. 3) reduces 
to a single term thereby yielding a useful recursion relation. 

Note that while the recursion relation is expressed in terms of the kinematic matrix S, 
the underlying Gram matrix G, also plays an important role. For example, for > 6 the 
Gram determinant vanishes and this has generic consequences for the parameters {yi}. To 
understand this we decompose the kinematic matrix. The equation to be solved is [16], 

N N (\ 1 \ 

SijVj = cx(1,1,...,1)e^ -Gavj + -Gjjyj - G ijyj J (A.4) 

3=1 3=1 V ' 

for both c = and c = 1. In other words, 

N-i j N-i N 

Yl Gi j y j = 2 Bwu ^2 WjVj = c ' Y y j = B ( A - 5 ) 

3=1 3=1 3=1 

where W{ = \Ga. 
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Using momentum conservation at the level of the loop integral, we can enforce B = 
det(G) = by setting y^ = —Yl^iVj- ^ n fact, this particular choice guarantees that a 
solution exists and that, 

N-l N-l N 

Yj = °' Yj w m = c ' b = Y y j = °- ( A - 6 ) 

j=l j=l j=l 

Therefore any vector in the kernel of G satisfying w.y = c is a solution of (S.y)i = c for 
N > 6. 

We can now construct the specific {yi} for all the cases specified in sec. 2.3. To do so 
we apply the Singular Value Decomposition (SVD) technique. Using the decomposition we 
can explicitly construct the parameters {yi} for Eqs. (2.18), (2.28) and (2.30). Note that 
Eq. (2.19) follows from Eq. (2.18) by summing over index k and applying identity (A. 2). 
Subsequently, Eq. (2.21) follows by combining Eqs. (2.19) and (2.18). 

The three distinct situations are: 

I. N < 6: det(S) + 0, £f=i S lj b j = 5 ik 

In this case the inverse of the kinematic matrix exists and the range of the matrix 
covers all of parameter space, specifically the unit vectors. 

While many algorithms exist for calculating the inverse, we still use the stable SVD 
for the kinematic matrix. This gives good control over the numerical accuracy of the 
matrix inversion which will be important as the external momenta approach exceptional 
configurations. The SVD is given by 

N 

Sij = Y^kUikVkj, (A.7) 
k=i 

where the are the non-zero eigenvalues of S and the matrices u and v are orthog- 

onal i.e. (v.v T )ij = (u.u T )ij = Sij. For an explicit algorithm which calculates both the 
eigenvalues and the orthogonal matrices see e.g. [34]. 
For any k < N, the choice, 

N 

bi = S- k 1 = Y"I lv v u Jk, ( A - 8 ) 
j'=i 

yields <%6j = <5jfc which immediately leads to Eq. (2.18). Note that because of Eq. (A.6), 
when N = 6 we find the special identity B = Ylij S^ 1 = 0, which then leads directly to 
Eq. (2.26). 

II. N > 7: det(S) = 0, £)JLi S ijZj = 

In this case the inverse is no longer defined. However for any vector in the kernel of 
the singular matrix S we have Ylj S-ij z j = 0- By selecting Zi with the property J2i z i = 0' 
we find Eq. (2.28). 
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To construct this solution, we simply solve Eq. (A. 6) by picking a vector out of the 
kernel of G such that w.z = For explicit construction of these vectors we first apply the 
SVD to the Gram matrix (and not the kinematic matrix), 



G v = ^Z^kUikvly (A.9) 



k=i 

With this decomposition, one of the possible choices for Zi is 



W e v i5 - W 5 v ie . 

Zl = w 5 + w 6 ' (A - 10) 



where 

N-l 



2 



Note that this solutions has the special property B = J2iLi z i = by construction. 
III. TV > 7: det(S) = 0, £V Sijrj = 1 

The fact that the vector (1, 1, . . . , 1) is in the range of the singular matrix is a highly 
non-trivial statement. It is due to the properties of the underlying decay/scattering process. 
With the property that ■ SijVj = 1 and Y2j r j = we readily obtain Eq. (2.30). 

By inspecting Eq. (A. 6), we know that a solution exists. We simply pick a vector out 
of the kernel of G such that w.r = 0. 

A explicit construction of the parameters r\ again passes through the SVD of the Gram 
matrix, Eq. (A.9). An example of the rj satisfying all of the above requirements is given 
by, 

"-w.- (A ' 12) 



With these explicit constructions of the solutions for (b,z,r), we have fully specified 
the recursion relations. The SVD yields, not only an explicit construction of the solution, 
but also a diagnostics tool for exceptional momentum configurations. These configurations 
can occur, even in non-singular cases. While integrable, they may require special treatment 
to maintain numerical accuracy. The SVD is particular useful to identify these situations. 

B. Analytic forms for the divergent integrals 

The one-loop bubble integral with external momentum scale Q\ with arbitrary powers of 
massless propagators is given by 

I 2 (D,{v u v 2 })-(-l) (-Q,), r(i/l)rh)r(I) _ ff) 
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Inserting specific choices for {vi, v 2 }, we see that, 
J 2 (D = 4-2e;{l,l}) = c r ' 



e(l - 2e) 



(-Qi) , 



7 2 (D = 6-2e;{l,l}) = c r 



I 2 (D = 6-2 e; {2,l}) = -c r 
/ 2 ( J D = 8-2e;{2,2}) = c r 
/ 2 (D = 8-2e; {3,1}) = c r 



2e(l-2e)(3-2e) 



(-0?) 



2e(l - 2e) 

(1-6) 



where 



cr 



2e(l - 2e)(3 - 2e) 
4e(l-2e)(3-2e) 
r(l-e) 2 r(l + e) 



(-<??)"'■ 



(B.2) 
(B.3) 
(B.4) 
(B.5) 
(B.6) 

(B.7) 



r(l-2e) 

The one- loop triangle integral with two on-shell legs, one off-shell external momentum 
scale Q\ and with arbitrary powers of massless propagators is given by 



h(D;{ Vl ,v 2 ,uz}) = {-iy{-Ql)T 



r(f-, 1 -, 2 )r(f-, 1 -, 3 )r(a-f) 



(B.8) 



r(iA2)r(iA,)r(£>-<7) 

Note that v\ refers to the propagator opposite the off-shell external leg. 

The one-loop triangle integral with one on-shell legs, two off-shell external momentum 
scales Q\ > Q 2 and with arbitrary powers of massless propagators is given by 



h{D;{v 1 ,v 2 ,vz}) 
= {-If (-Q\) k > 



(B.9) 



2 ,oT (f - vi ~ vi) r(f -vi- ^)Y{p - f ) 



i> 2 )r> 3 )r(D - a) 

D 

x 2 Fi [u!, a - — , 1 + 1^1 + 1/3 



Of 
2'Q? 



r(i/i)r(^)r(D - a) 

Note that ^3 refers to the propagator opposite the on-shell external leg. 



D 

x 2F1 ( i/ 2 , — - 1/3, 1 + -5- - - ^3, 7^2 



(B.10) 



C. The scalar iV-point function 

Using Eq. (3.15) can write down the singular structure for the scalar TV-point function with 
all legs on-shell, 



C(12-JV) 



N 



3=3 



12J 
N 



(C.l) 
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Integrating out the loop momentum in the scalar triangle integrals yields 




N 



(823-j ~ S 3-j) 




I\k=3,k^j( S 23-j s 3-k - S 3 ...jS23-k) 



((-S23- 3 r ~ (SS-iD ■ (C2) 



This expression agrees with the singular structure singularity structure for the massless six- 
point function with all momenta on-shell given in Ref. [16] and the corresponding massless 
five-point function with all momenta on-shell given in Ref. [32] as well as the massless 
four-point function of Ref. [14] . 

It is straightforward to use Eq. (3.15) to check the known singularity structure for 
massless box graphs with all possible combinations of off-shell conditions for the external 
momenta given in Ref. [14] as well as the massless pentagon with one off-shell external 
momenta [32]. 
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